PSEUDO-TRANSIENT  CONTINUATION  FOR  NONSMOOTH  NONLINEAR 

EQUATIONS  * 
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Abstract.  Pseudo-transient  continuation  is  a  Newton-like  iterative  method  for  computing  steady-state  solutions 
of  differential  equations  in  cases  where  the  initial  data  is  far  from  a  steady  state.  The  iteration  mimics  a  temporal 
integration  scheme,  with  the  time  step  being  increased  as  steady  state  is  approached.  The  iteration  is  an  inexact 
Newton  iteration  in  the  terminal  phase. 

In  this  paper  we  show  how  steady-state  solutions  to  certain  ordinary  and  differential  algebraic  equations  with 
nonsmooth  dynamics  can  be  computed  with  the  method  of  pseudo-transient  continuation.  An  example  of  such  a  case 
is  a  discretized  partial  differential  equation  with  a  Lipschitz  continuous,  but  non-differentiable,  constitutive  relation  as 
part  of  the  nonlinearity.  In  this  case  we  can  approximate  a  generalized  derivative  with  a  difference  quotient. 

The  existing  theory  for  pseudo-transient  continuation  requires  Lipschitz  continuity  of  the  Jacobian.  Newton-like 
methods  for  nonsmooth  equations  have  been  globalized  by  trust-region  methods,  smooth  approximations,  and  splitting 
methods  in  the  past,  but  these  approaches  require  problem-specific  components  in  an  algorithm.  The  method  in  this 
paper  addresses  the  nonsmoothness  directly. 
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1.  Introduction.  In  this  paper  we  show  how  pseudo-transient  continuation  ('Ttc  )  can  be  used 
to  solve  a  class  of  nonsmooth  nonlinear  equations,  'btc  is  a  predictor-corrector  method  for  efficient 
integration  of  a  time-dependent  differential  equation  to  steady  state.  The  objective  of  the  method  is 
not  temporal  accuracy,  but  rather  to  resolve  the  transient  behavior  of  the  solution  until  the  iteration 
is  close  to  steady  state,  and  then  to  increase  the  “time  step”  and  transition  to  a  fast  Newton-like 
method. 

In  this  paper  we  extend  the  theoretical  convergence  results  of  [7, 18]  to  problems  with  cer¬ 
tain  nonsmooth  nonlinearities  and,  thereby,  partially  explain  the  results  reported  in  [8, 10].  We  also 
show  how  generalized  derivatives  can  be  approximated  by  finite  differences,  and  how  those  approx¬ 
imate  derivatives  can  be  used  effectively  both  in  locally  convergent  iterations,  such  as  those  which 
arise  in  temporal  integration,  and  in  the  context  of  4Tc  .  This  aspect  of  the  work  is  motivated  by 
several  papers  on  simulation  of  unsaturated  flow,  [10, 14, 15,24,30,31],  in  which  Lipschitz  contin¬ 
uous  spline  approximations  to  the  non-Lipschitz  continuous  van  Geneuchten  and  Mualem  [25,33] 
constitutive  laws  are  used.  These  nonsmooth  functions  are  then  differentiated  with  finite  differ¬ 
ences  as  if  they  were  smooth.  The  results  in  this  paper  explain  the  success  reported  in  those 
papers.  Another  aspect  of  the  paper  is  an  extension  of  the  local  results  in  [9, 21, 27, 28]. 

'Ttc  methods  are  particularly  appropriate  for  the  types  of  nonsmooth  nonlinearities  which  we 
discuss  in  this  paper.  Traditional  methods  for  globalizing  iterative  methods  for  nonlinear  equa- 
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tions,  such  as  line  searches,  can  fail  as  commonly  implemented  in  practice  for  both  smooth  and 
nonsmooth  equations  [7, 8, 18].  The  existing  global  convergence  results  for  nonsmooth  nonlinear 
equations  are  either  based  on  line  searches  for  a  inexact  Newton  formulation  [9,22],  a  sequence 
of  smooth  approximations  [4,29],  or  explicit  treatment  of  the  nonsmoothness  [20].  Only  the  latter 
admits  approximation  of  the  generalized  Jacobian  by  a  difference,  'Ttc  allows  one  to  deal  with  the 
nonsmoothness  directly  and  exploits  the  dynamics  to  guarantee  convergence  to  x*,  the  solution  one 
wants. 

In  the  remainder  of  this  introductory  section  we  review  the  relevant  results  from  nonsmooth 
analysis  (§  2.1)  and  Ttc  (§  2.2).  Then  we  describe  the  setting  for  the  new  results.  In  §  3  we 
show  how  finite  difference  approximations  of  generalized  Jacobians  affect  the  local  convergence 
of  inexact  Newton  methods  for  nonsmooth  problems.  We  use  those  results  in  §  4,  where  we  state 
and  prove  our  local  and  global  convergence  results  for  Ttc  .  We  present  a  numerical  example  in 
§5. 

Some  extensions  of  our  results  to  infinite  dimensions  are  possible,  using  ideas  from  [5, 11, 
12, 19, 32]  if  the  appropriate  compactness  conditions  hold.  These  extensions  will  be  explored  in  a 
subsequent  paper. 

2.  Previous  results.  In  this  section  we  review  the  prior  results  about  Ttc  and  nonsmooth 
analysis  that  we  will  need  for  this  paper. 

In  this  paper  the  norm  will  be  the  a  scaled  discrete  l2  norm  on  RN, 


unless  stated  otherwise.  The  ball  of  radius  e  about  a  point  x  G  RN  will  be  denoted 

B(x,  e)  =  {z  |  ||a;  —  z\\  <  e}. 

As  is  standard,  we  will  let  x*  denote  solution  of  F(x )  =  0,  and 

e  =  x  —  x* 


the  error.  We  will  let  (x)i  denote  the  / 1 h  component  of  the  vector  x. 

2.1.  Nonsmooth  Analysis.  In  this  section  we  review  the  concepts  from  nonsmooth  analysis 
[6,23]  which  we  will  need  for  our  convergence  results.  We  then  state  the  local  convergence  result 
from  [21,28]  which  extend  in  this  paper. 

Let  F  :  RN  — >  RN  be  locally  Lipschitz  continuous.  This  implies  that  F  is  Frechet  differen¬ 
tiable  almost  everywhere,  and  that  the  directional  derivatives 


F'(x  :  w)  —  lim 

v  ’  h^O 


F(x  +  hw)  —  F(x) 


exist  for  all  x,  w  G  RN . 

We  let  DF  denote  the  set  of  points  where  F  is  Frechet  differentiable.  The  generalized  Jacobian 
[6]  of  F  at  u  G  RN  is  the  set 
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where  co  denotes  the  convex  hull. 

We  will  consider  extensions  of  the  Newton-like  iteration 


(2.2)  xn+i  =  xn  -  Vn  1F(xn) 

where,  Vn  G  d F(xn),  and,  as  is  standard  xn  is  the  current  approximation  to  a  solution  x*  and  xn+  \ 
the  new  approximation.  Our  results  will  be  stated  in  terms  of  an  inexact  formulation, 

(2.3)  Xn+1  =  xn  +  s, 


where 

(2.4)  \\Vns  +  F(xn)\\  <  77n||F(xn)||, 
and  Vn  G  dF(xn). 

The  most  important  concept  is  that  of  semismoothness  [23,28]. 

Definition  2.1.  F  is  semismooth  at  x  G  RN  if  F  is  locally  Lipschitz  continuous  and  for  all 
w  G  RN,  the  limit 

(2.5)  lirn  {Vw'\ 

V  EdF(x-\-tw'),w'^w,tiO 

exists. 

Semismoothness  is  a  useful  concept  [5,28,32]  in  the  proofs  of  convergence  and  local  conver¬ 
gence  rates  of  the  iteration  (2.2).  In  the  standard  theory  for  Lipschitz  continuously  differentiable  F, 
local  quadratic  convergence  follows  from  nonsingularity  of  the  Jacobian  F'(x*)  at  the  solution  and 
the  Lipschitz  continuity  of  F'.  In  the  nonsmooth  case,  one  must  prove  that  the  Newton  iteration  is 
well  defined  and  quantify  the  degree  of  nonsmoothness  to  obtain  convergence  rates. 

Lemma  2.2,  taken  from  [27],  and  Lemma  2.4,  are  the  results  that  are  needed  to  prove  local 
superlinear  convergence  of  (2.2). 

Lemma  2.2.  F  is  semismooth  at  x  G  RN  if  and  only  if 


(2.6) 


lim 

w — >0,  V  £dF(x+w) 


||  F(x  +  w)  —  F(x)  —  Vw  || 

M  " 


o. 


To  obtain  convergence  rates  one  needs  a  stronger  condition  than  semismoothness  [28]. 

Definition  2.3.  F  is  semismooth  of  order  p  at  x  if  for  all  w  G  RN  and  V  G  d  F(x  +  w) 


(2.7) 


F(x  +  w)  —  F(x)  —  Vw  =  0(||w||1+p) 


as  w  — »  0. 

LEMMA  2.4.  Let  F  be  semismooth,  F(x*)  =  0,  and  assume  that  all  matrices  in  dF(x*) 
are  nonsingular.  Then  there  are  M  and  A  such  that  if  x  G  B(x*,  A)  and  V  G  dF(x),  then 
||  V-1 1|  <  M. 

These  results  have  been  used  to  prove  several  convergence  theorems  [9, 21, 27, 28]  for  (2.2) 
and  (2.3).  Theorem  2.5  is  a  combination  of  the  local  convergence  results,  and  is  the  basis  for  the 
new  algorithms  in  this  paper. 

Theorem  2.5.  Let  F  :  RN  — >  RN  with  F(x*)  =  0.  Assume  that  F  is  semismooth  at  x*. 
Then  there  are  fj,  S.  K  >  0  such  that  ifx0  G  B(x* .  5)  and  tjn  <  fj  then  the  inexact  Newton  iteration 
(2.3)  converges  to  x*  and 

H^n-i-ill  A  -Ar/n 1 1 en 1 1  T  o(||^n||)- 
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Moreover,  if  F  is  semismooth  of  order  p  at  x*,  then 


||en+i||  <  K(rjn\\en\\  +  ||en||1+p). 


In  previous  work  [11,16, 19,20]  on  nonsmooth  nonlinear  equations  in  function  spaces  and  their 
discretizations,  we  used  properties  of  the  solution  to  isolate  a  smooth  component  of  the  nonlinear¬ 
ity.  Each  problem  required  a  different  approach,  and  all  assumed  that  the  nonsmooth  component 
was  small.  In  those  papers,  mesh-independent  convergence  results  were  obtained,  and  standard 
implementations  of  matrix-free  Newton-Krylov  methods  worked  well. 

The  formulation  we  consider  in  this  paper  is  different.  Here  one  does  not  have  to  explicitly 
split  the  operator  into  smooth  and  nonsmooth  parts,  a  significant  advantage  in  complicated  appli¬ 
cations  [30].  However,  we  know  of  no  general  proofs  of  mesh-independent  convergence  rates,  a 
problem  also  mentioned  in  [32].  In  fact,  the  numerical  results  in  §  5  show  mesh-dependent  perfor¬ 
mance  of  the  iteration,  especially  in  the  mid-range.  Mesh-dependent  convergence  was  also  reported 
in  [4].  In  §  5.3  we  illustrate  this  phenomenon  and  show  how  a  nested  iteration  can  overcome  it. 

Numerical  differentiation,  a  topic  we  consider  in  §  3.1,  is  a  simple  matter  if  the  smooth  and 
nonsmooth  parts  of  the  nonlinearity  are  split.  Here,  we  can  prove  accuracy  only  if  one  is  differenti¬ 
ating  in  coordinate  directions,  and  only  then  for  special  classes  of  operators.  Since  the  directions  in 
a  matrix-free  Newton-Krylov  solver  are  not  predictable,  our  results  do  not  apply  to  those  methods. 

2.2.  Pseudo- Transient  Continuation.  The  objective  of  Ttc  ,  as  we  present  it  here,  is  to  find 
the  steady  state  solution  of  the  semi-explicit  index-one  differential  algebraic  equation  (DAE) 

<2.8)  °( “)'= -($:»> )=-PW>  x(0)=xo- 

Here  x  =  ( u,v)T  G  C([0,  oo],  RNl+N*.  The  functions  u  :  [0,  oo]  — >  RNl  and  v  :  [0,  oo]  — >  RN2 

are  to  be  found.  The  differential  variables  u  and  the  algebraic  variables  v  are  clearly  separated  in 
the  semi-explicit  case  where 


where  Du  is  a  nonsingular  scaling  matrix.  A  good  general  reference  for  DAEs  is  [3]. 

We  assume  that  the  initial  data  for  (2.8)  are  consistent  (;.  e.  g(u(0),v(0))  =  0)  and  seek  the 
solution  x*  to  F(x*)  =  0  that  satisfies 


lim  x(t)  =  x*. 

t— XX) 

If  (2.8)  is  a  discretization  in  space  of  a  PDE,  and  the  initial  data  is  far  from  the  desired  steady 
state,  the  application  of  a  conventional  method,  such  as  a  line  search  [17],  to  the  time-independent 
equation 

F(x)  =  0, 

may  fail  to  converge.  Possible  failure  modes  [8]  are  stagnation  of  the  iteration  at  a  singularity  of 
F',  the  Jacobian  of  F,  or  or  finding  a  solution  other  than  x*. 
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We  formulate  Ttc  as 

(2.9)  xn+i  =xn-  +  F\xn))~1F(xn). 

In  (2.9),  {<5n}  is  adjusted  to  efficiently  find  the  steady  state  solution  rather  than  to  enforce  temporal 
accuracy. 

The  convergence  results  in  [7, 18]  assume  that  the  time  step  was  updated  with  “switched  evo¬ 
lution  relaxation”  (SER)  [26],  i.  e. 

,  . (s  lin*»-l)ll  r  'l 

(2.10)  0n  max  I  0n—i  .I  ,  x  ||  i^max  I  • 

In  [7]  we  proved  convergence  for  smooth  F  under  the  assumptions  that  the  DAE  has  index 
one  in  a  certain  uniform  sense,  has  a  global  solution  in  time,  and  that  the  solution  converges  to  a 
steady  state.  The  convergence  result  for  the  exact  rj  —  0  case  is 

(2.11)  \\xn+1  -  x*\\  =  0( \\xn  -  a* || (S~lx  +  \\xn  -  a:*||)) 
as  n  — >  oo. 

In  this  paper  we  relax  the  smoothness  assumptions  and  consider  the  iteration 

(2.12)  xn+i  =xn  +  s, 
where 

(2.13)  IKS-'D  +  V(xn))s  +  F(xn)\\  <  Vn\\F(xn)\\ 
where  V (xn)  is  near  to  the  set  d F(xn)  the  sense  that 

(2.14)  V (xn)  G  T>(xn,  C,  h),  for  some  small  h. 
where 

(2.15)  V{xn,C,h)  =  {V\\\V  -V\\  <  Ch,  for  some  V  e  dF{x)  and  y||  <  Ch}. 

The  sense  in  which  V  (xn)  is  close  to  dF(xn)  is  technical  because  dF{x)  is  not  continuous  in 
x.  The  requirement  that  V ( xn )  e  F>(xn,  C,  h )  is,  in  a  sense,  a  requirement  that  a  combination  of 
the  forward  and  backward  error  be  small. 

In  [7,  8, 10, 18]  we  report  on  compuational  results  that  show  that  both  the  local  and  global 
phases  of  the  'Etc  iteration  perform  as  (2. 1 1)  predicts  even  if  the  nonlinearity  is  nonsmooth  [7,8, 10] 
and  the  derivative  is  approximated  by  differencing  [10, 13, 14,30]. 

3.  Local  Convergence  Theory.  In  this  section  we  analyze  the  accuracy  of  a  finite  difference 
approximation  of  a  generalized  Jacobian  in  the  case  where  the  nonsmoothness  arises  from  a  sub¬ 
stitution  operator.  The  approximation  is  accurate  in  a  combined  forward  and  backward  sense,  and 
this  affects  not  only  the  convergence  speed  of  an  inexact  Newton  iteration,  but  also  the  limiting 
accuracy. 
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3.1.  Finite  difference  approximations.  The  results  in  this  paper  are  motivated  in  part  by  our 
experience  with  nonsmooth  nonlinear  substitution  operators.  A  substitution  operator  on  R A  has 
the  form 

(3.1)  $(x)  =  (0(zi),  ■  ■  ■  ,f(xN))T 

where  0  :  R  — >  R.  The  generalized  Jacobian  of  <f>  is  the  set  of  diagonal  matrices 

d$(x)  =  (<90(zi), . . . ,  df(xN))T. 

In  this  section,  we  consider  maps  that  are  compositions  of  smooth  maps  with  substitution 
operators.  Let 

G(x)  =  S($(x)) 

where  S  is  differentiable.  The  smoothness  of  S  and  the  definition  (2.1)  of  d G  imply  that 

d  G(x)  =  S'mx))d${x) 

where  S'  is  the  Jacobian  of  S.  Of  interest  here  is  the  approximation  of  ()G  with  a  finite  difference 
approximation  using  the  coordinate  directions. 

Let  d£G{  x)  be  the  matrix  whose  ?’th  column  is 

G(x  +  hli)  -  G(x ) 
h  ’ 

where  1,:  is  the  unit  vector  in  the  /'th  coordinate  direction.  We  show  that  the  forward  difference 
df^F  approximates  d G(x)  in  the  sense  described  by  (2.15). 

Theorem  3.1.  Let  0  :  R  — >  R  be  Lipschitz  continuous  and  differentiable  except  at  finitely 
many  points  Let  S  be  Lipschitz  continuously  differentiable  in  RN .  Then  there  is  C  >  0 

such  that  for  all  h  sufficiently  small 

(3.2)  d %G{x)  eV(x,C,h). 


Proof.  We  begin  by  showing  that  it  suffices  to  prove  the  result  for  scalar  functions.  Differen¬ 
tiability  of  S  and  the  Lipschitz  continuity  of  0  imply  that 

G(x  +  hli)  -  G(x)  =  S’mx))($(x  +  hh)  -  $(x))  +  0(h2) 

=  S"($(x))($(x  +  hit)  -  $(x))  +  0(h2) 

for  all  x  such  that  ||a;  —  x \  <  h.  Hence  we  need  only  prove  the  result  for  substitution  operators. 
Since  <f>  is  a  substitution  operator,  the  ith  component  of  <f>(a;  +  hli)  —  Q(x)  is 

0((®)i  +  h)  -  0((x)i), 

and  we  need  only  consider  scalar  functions.  Now  let, 


h  <  min  ||^  -  ^Hoo, 
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then  at  most  one  £  is  in  the  interval  |  (x)t,  (x)l-\-h\.  If  o  is  differentiable  in  the  interval  |  (x)t,  (x)i+h] 
then 


4>{{x)i  +  h)  -  4>{{x)i) 
h 


0'((aO*)  +  0(h) 


and  we  let  the  7th  component  of  x  be  (x)*. 

Now  assume  that  £,•  G  [(x)j,  (x),  +  /i]  for  some  j.  Let  <j>'+(£j)  and  0'_  (£.,■)  be  the  right  and  left 
handed  derivatives  at  £j 


0±&) 


lim  0(£j  ±  h)  -  0(£j)) 
h— >0  ±/i 


Let  -  (x)j  =  i//t,  for  v  G  [0,1].  Since 


0((x)i  +  h)  4>((x)i)  =  0(£j  +  (1  -  u)h)  -  0(£j  -  vh) 


=  4>(£j  +  (1  “  v)h)  -  0(0)  +  0(0)  -  0(0  -  uh) 

=  (1  -  ^)0'+(O)  +  ^0-(O)  +  0(h2). 

Since 

(1  -  Z/)0+( 0)  +  ^0-(O)  e  ^0(0) 

for  all  z/  G  [0, 1],  the  proof  is  complete  with  x%  =  □ 

A  similar  result  holds  for  central  differences.  Let  df-'G(x)  be  the  matrix  whose  7 1 h  column  is 

G(x  +  hli)  —  G(x  —  hli) 

2 h  ‘ 

If  S  is  Lipschitz  continuously  twice  differentiable  and  (j)  is  piecewise  Lipschitz  continuously  twice 
differentiable,  then  the  statement  of  Theorem  3.1  with 

||I/-y||  <  Ch 

in  (2.15)  replaced  by 

(3.3)  \\V-V\\<Ch2. 

3.2.  Local  Convergence.  If  the  generalized  Jacobian  is  approximated  by  a  finite  differece,  one 
cannot  expect  asymptotic  convergence,  because  the  accuracy  in  the  terminal  phase  of  the  iteration 
will  be  limited  by  the  accuracy  in  the  derivative.  We  quantify  this  in  Theorem  3.2,  which  extends 
the  existing  local  convergence  theorems  for  inexact  Newton  methods  for  semi- smooth  equations. 
The  new  assumption  that  V (x)  G  V(x,  C,  h )  is  motivated  by  the  results  in  §  3.1. 

THEOREM  3.2.  Assume  that  F  is  semismooth  at  x*,  F(x*)  =  0,  and  that  all  matrices  in 
d F(x*)  are  nonsingular.  Assume  that  there  is  C  >  0  such  that 

(3.4)  V(x)  G  V(x,C,h) 
for  all  x  sufficiently  near  x*. 

Then  there  is  5  such  that  ifx o  G  B(x*,  5),  {t]n}  and  h  are  sufficiently  small,  then  the  iteration 


(3.5) 


Xn+l  Xn  T  S , 
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where 

(3.6)  \\V(xn)s  +  F(xn)\\<V\\F(xn)l 
converges  to  x*.  Moreover,  there  is  K  >  0  such  that 

(3.7)  ||e™+i||  <  K((Vn  + h)\\en\\  +  h)  +  o(\\en\\)7 
or,  if  F  is  is  semismooth  of  order  p  at  x*,  then 

(3.8)  ||en+i||  <  K((r)n  +  h)\\en\\  +  ||en||1+p  +  h). 

Proof.  The  plan  of  the  proof  is  to  compare  xn+\  with  the  Newton  iteration  from  xn,  where  xn 
is  the  point  specified  in  the  definition  of  V.  We  can  then  apply  Theorem  2.5. 

Let  5  and  h  be  small  enough  so  that 

(3.9)  ||f1||  <  M,  for  all  u  e  B{u* ,  h  +  5), 

which  we  can  do  by  Lemma  2.4.  We  assume  that  xn  6  B(x *,  5)  and  will,  reducing  S  and  h  if 
necessary,  show  that  xn+\  G  B(x*,  5)  and  that 

By  assumption,  there  are  xn  6  B(xn,  h )  and  Vn  e  d F(xn)  such  that 

\\V(xn)-Vn\\<Ch. 

Hence  the  step  s  is  nearly  an  inexact  Newton  step  from  xn. 

By  (3.9),  for  h  sufficiently  small, 

||^«T1||  <  1  /(AT1  -  Ch)  <  2 M 

and  hence 

||s||  <  2 M(rjn  +  l)||F(a;n)||. 

Therefore, 

||Ks  +  F(xn) ||  <  || V(xn)s  +  F(xn) ||  +  Ch\\s\\ 

(3.10)  <  ||y(xn)s  +  F(xn)  ||  +  \\F(xn)  -  F(xn)  ||  +  Ch\\s\\ 

<  r]n\\F{xn)\\  +  Lh  +  Ch(2M{\  +  rin)\\F(xn)\\, 

where  L  is  the  Lipschitz  constant  of  F.  Since 

\\F(xn)\\<\\F(xn)\\+Lh, 

we  may  set 

Kq  =  1  +  2 L  +  2MC(l  +  r]n)  <1  +  2 L  +  AM C, 

and  obtain 

(3.11)  \\Vns  +  F{xn)\\  <K0((rin  +  h)\\F(xn)\\  +  h). 
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The  inexact  Newton  condition  (3.6)  and  (3.11)  imply  that 

xn+l  =  xn  +  V~1(F(xn)  +  r  n) 

with 

||r„||  <  Kobrin  +  h)\\F(xn)\\  +  h ). 

Hence, 

||e„+i||  =  \\en  +  V~lF(xn)  +  14_1rn|| 

<  MK0((r)n  +  /t)||F(xn)||  +  h)  +  o(||en||). 

If  F  is  semismooth  of  order  p  at  x*.  Theorem  2.5  implies  that  there  is  K1  >  0  such  that 

IK+ill  =  \\en  +  V-'Fixn)  +  V~lrn\\ 

<  ^i||e„H1+p  +  M K0{j]n\\F (xn)\\  +  h). 

Since  ||F(xn)||  <  L||en||  and  ||en||  <  ||en||  +  h,  we  obtain  (3.8)  with 

K  =  2/C,  +  MK0(  1  +  L), 


and  complete  the  proof.  □ 

3.3.  Optimal  choice  of  h.  If  xn  7^  xn,  then  the  estimate  (3.7)  and  (3.8)  do  not  imply  con¬ 
vergence,  but  stagnation  once  the  error  is  0(h).  This  is  analogous  to  convergence  results  [17]  for 
Newton’s  method  when  there  are  errors,  such  as  floating  point  roundoff,  in  the  evaluation  of  F.  In 
this  case,  however,  h  is  larger  than  floating  point  roundoff,  and  we  can  combine  Theorems  3.1  and 
3.2  to  estimate  the  optimal  choice  of  h. 

Suppose  that  F  is  piecewise  C 1  (and  hence  semismooth  of  order  1  [23])  and  can  be  evaluated 
up  to  an  absolute  error  of  €f-  If  we  incorporate  the  error  in  F  into  the  result  of  Theorem  3.1  in  the 
standard  way  [17]  we  obtain 

d%G(x)  eV(x,C',ti) 

where  h!  =  0(h  +  ep/h).  Then  the  estimate  (3.8)  becomes 

(3.12)  ||en+i||  ^  K ((Vn  +  h  +  e^//i)||en||  +  ||en||2  +  h). 

If  we  solve  the  equation  for  the  step  exactly,  then  r/n  =  0.  In  that  case,  if  ||e„||  =  0(h^2), 

then 

(3.13)  Hen+l||  =  o  • 

The  term  on  the  right  of  (3.13)  is  minimized  when 

(3.14)  h  =  0(ep3). 

If,  for  example,  ep  ~  10~15  is  double  precision  floating  point  roundoff,  (3.14)  would  say  that  the 
best  results  would  be  obtained  if  h  &  10-10,  rather  that  10~8  as  a  conventional  analysis  would 
predict.  We  provide  numerical  evidence  for  this  in  §  5.2. 
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4.  Convergence  of  TTc  .  The  analysis  of  Ttc  in  this  paper  follows  the  pattern  of  [7,  18], 
considering  the  iteration  in  two  phases.  For  phase  one,  the  initial  or  global  phase,  we  show  that 
Ttc  is  a  consistent  convergent  scheme  for  integration  of  the  DAE.  The  scheme  will  be  first  order 
if  F  is  semismooth  of  order  1,  order  1  +  p  if  F  is  semismooth  of  order  p  <  1,  and  convergent,  but 
with  no  order,  if  F  is  merely  semismooth. 

From  the  analysis  of  the  global  phase  we  will  conclude  that,  for  sufficiently  small  <50,  the 
iteration  will  approach  x*.  For  the  second  local  phase  of  the  iteration,  we  show  that  if  x  is  near  x* 
and  {<5n}  is  bounded  away  from  zero,  then  8n  — >  Smax  and  then  the  terminal  phase  of  convergence 
can  be  described  by  the  results  in  §  3. 

The  analysis  of  the  local  phase  does  not  depend  on  the  dynamics,  and  we  will  defer  the  detailed 
assumptions  on  the  DAE  until  §  4.2. 

4.1.  Local  Phase.  We  consider  the  local  phase  first,  as  we  did  in  [7, 18],  in  order  to  establish 
targets  for  the  integration  in  the  global  phase.  We  seek  to  find  eL  so  that  if  x0  G  B(x* ,  ef)  and  {()„  } 
is  bounded  away  from  zero,  then  {a;n}  and  {5n}  in  (2.12)  satisfy  Xn  *  X  and  8n  *  8max  or 

The  local  convergence  rates  in  the  terminal  phase  depend  on 

ASSUMPTION  4.1.  F  is  semismooth  at  x*.  There  are  C.  h,  /3,€l  >  0  such  that  for  all  x  G 
B(x*,  €l)  and  all  8  >  0 

\\(D  +  8V(x))-1D\\<l/(l  +  ^8), 

and 

V{x)  G  V(x,C,h). 


THEOREM  4.1.  Let  the  assumptions  of  Theorem  3.2  and  Assumption  4.1  hold.  Let  { 8n }  he 
given  by  (2.10).  Then  there  are  Ct  and  eT  such  that  if  {rjn}  is  sufficiently  small  and  x$  G  B(x* ,  €t), 
then  either  infn  8n  =  0  or  8n  — >  8max,  the  yVtc  iteration  converges,  and,  for  n  sufficiently  large 

(4.1)  ||en+i||  <  CT{fj]n  +  8n  1  +  h) ||en||  +  h)  +  o(||en||) 
or,  if  F  is  semismooth  of  order  p, 

(4.2)  ||en+i||  <  CT(||en||1+:p  +  ft]n  +  8nl  +  /i)||en||  +  IT). 


Proof.  We  assume  that  a;0  is  near  enough  to  x*  so  that  the  conclusions  of  Theorem  3.2  hold.  If 
xn  G  BUt),  then  following  the  proof  of  Theorem  3.2, 

en+ i  =  Cn  ~  (8n  1 D  +  Vn)  1F(xn)  +  rn 


where 

1 1  t'n  1 1  =  0((rin  +  h)\\F(xn)\\  +  h). 

Semi  smoothness  and  our  assumptions  imply  that 

F{xn)  -  Vnen  =  0(h)  +  o(||en||) 
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and  hence 

c-n+ 1  cn  (3n  D  +  Vn)  Vnen  -\-  Rn 
=  (S^D  +  Vn)~1S~1Den  +  Rn, 

where 

Rn  =  0((rjn  +  h)\\F(xn)\\  +h)  +  o(||en||). 

If  5n  >  5*  for  all  n,  then  Assumption  4.1  implies  that 

\\(5-lD  +  Vn)-H-1D\\<l/{l  +  (35*). 

This  implies  that  the  iteration  is  q-linearly  convergent,  and  hence  8n  — >  Smax  and  xn  — >  x*. 

The  completion  of  the  proof  for  large  5n  is  a  direct  consequence  of  the  Theorem  3.2,  since  the 
inexact  Newton  condition 


||(iTB  +  V(i„))S  +  F(xn)||  <  7k||F(z„)|| 

implies  that  there  is  Ch  such  that 

||F(xn)s  +  F(xn)||  <  (vn  +  Chh)\\F(xn)\\  +  A.”1  ll-Dn || ||s|| , 

and  then  C  and  //,,  in  (3.10)  can  be  replaced  by  C  +  |  D\  \  |  and  //,,  +  8~ 1  +  CRi.  This  implies 
convergence  if  5max  is  sufficiently  large.  □ 

4.2.  Global  Phase.  In  the  analysis  of  the  global  phase  we  must  assume  that  the  4Tc  iteration 
is,  for  small  5,  a  stable  explicit  method  for  the  DAE  (2.8).  To  do  this  we  must  assume  that  the 
DAE  is  consistent  and  has  index  one.  In  the  smooth  case,  one  can  express  this  in  terms  of  the 
nonsingularity  of  gv,  the  Jacobian  of  g  with  respect  to  the  algebraic  variables.  In  the  nonsmooth 
case,  however,  one  must  take  the  limit  in  (2.1)  in  all  components  together.  This  means  that  the 
index  assumption  is  more  technical,  using  the  nonsingularity  of  the  matrix  pencil  d  ~ 1 D  +  V  ( x )  in 
part  5  of  Assumption  4.2. 

We  assume  that  V(x)  G  T>(x,C,h),  for  a  sufficiently  small  h.  We  decompose  operators 
V  G  OF  into  blocks 

(4.3)  V(x)  =  (  !)“  K”  )  , 

\  V  Vll  V  vv  J 

where  Vuu  G  duf, . . . ,  Vvv  G  dvg. 

With  this  in  mind  we  can  formulate  our  assumptions  on  the  dynamics.  Define  a  neighborhood 
of  the  trajectory  from  :r0 

(4.4)  S(e)  =  { z  |  inf  \\z  -  x(t) ||  <  e}. 

Assumption  4.2.  g(u0,  v0)  =  0,  i.  e.  the  initial  values  (tt0,  no)  are  consistent. 

There  are  Cq  G  (0,  67-/2),  where  6t  is  the  radius  from  Theorem  4.1,  such  that 

1.  F  is  semismooth  in  S(€q). 

2.  Vvv(x)  is  nonsingular  for  all  x  G  S(ea),  and  there  is  My  such  that  ||E),„(a;)_1||  <  My  for 
all  x  G  S(es). 
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3.  For  all  z0  G  S(es),  the  solution  of  Dz'  =  —F(z),  z(0)  =  z0  exists,  z(t)  G  S(ec)  for  all 

t,  and  lim^oo  z(t)  =  x*. 

4.  V  ( x )  G  T>(x ,  C,  h)  for  all  x  G  S(ec)- 

5.  Moreover,  there  are  Mo,  Mj  >  0  such  that  for  all  5  >  0, 

(a)  (5_1D  +  V  (x))  is  nonsingular  for  all  x  G  S(cg), 

(b)  ||(5_1.D  +  V(a:))||  <  MDfor  all  x  G  S(cg),  and 

(c)  ||  (d^D  +  ^(x))-1!!  <  Mo for  all  xG  S{eG). 

We  analyze  the  global  phase  by  simply  showing  that  the  global  truncation  error  of  the  scheme 

Xn+1  tCn  (5n  D  V  (xn))  f(ln)) 

is  of  order  p,  ie 

||^n  -  X{tn)  ||  =  0(5P) 

where  <5  =  max o<m<n5n.  This  will  imply,  similarly  to  [7, 18],  that  the  Ttc  iteration  will  correctly 
track  the  solution  until  xn  is  in  the  ball  of  local  convergence  required  by  Theorem  4.1. 

We  will  use  a  simple  consequence  of  pth  order  semi  smoothness. 

Lemma  4.2.  Let  Assumption  4.2  hold.  Let  x(t)  be  the  solution  to  (2.8).  Let  5  >  0  and  let 

^  =  x(t  +  6)  —  xit). 


(  <yu 
V  (Tv 


Then,  for  5,  h  sufficiently  small, 

(4.5)  (< 5~rD  +  V(x(t)))a  =  —F(x(t))  +  0(h)  +  o(S), 
and  if  F  is  semismooth  of  order  p, 

(4.6)  (d^D  +  V(x(t)))a  =  —F(x(t))  +  0(61+p  +  h), 
uniformly  in  t. 

Proof.  In  the  interests  of  brevity,  We  will  give  the  proof  for  h  —  0  and  F  semismooth  of  order 
p.  The  analysis  for  h  >  0  and  semismooth  F  follows  the  outlines  of  the  proofs  of  Theorems  3.2 
and  4.1. 

Write  x(t)  =  (u(t),  v(t))r.  By  integrating  the  DAE  (2.8),  we  see  that  v  is  a  Lipschitz  contin¬ 
uous  function  of  t.  The  nonsingularity  of  Vuu  implies  that  v  is  also  a  Lipschitz  continuous  function 
of  t.  This  Lipschitz  continuity  implies  that 


INI  =  0(5). 

Integrate  (2.8)  over  the  interval  [t,  t  +  5}  and  use  the  Lipschitz  continuity  of  F  to  obtain 

pt-\-S 

(4.7)  Da  =  —  J  F(x(r))dr  =  -5F(x(t  +  5))  +  0(52), 


uniformly  in  t. 
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By  the  definition  of  semismoothness  with  x  =  x(t  +  S)  and  w  =  —a  we  have,  for  5  sufficiently 
small, 

F(x(t  +  5))  =  F(x{t))  +  V(x(t))a  +  0(||tx|| 1+^) 

(4.8) 

=  F(x(t))  +  V  (x(t))a  +  0(51+p). 

The  estimate  (4.8)  is  uniform  in  t  because  the  set  {x(t)  \  t  >  0}  is  compact. 

Hence,  multiplying  (4.8)  by  delta  and  substituting  into  (4.7), 

(D  +  SV  (x(t)))a  =  -SF(x(t))  +  0(S2+P) 

and  the  proof  is  complete.  □ 

Lemma  4.2  will  imply  convergence  of  Ttc  in  the  same  way  as  in  the  smooth  case  [7].  The 
objective  is  to  show  that  for  <50  sufficiently  small,  the  Ttc  iteration  remains  in  the  tube  S(cc)  We 
state  the  result  and  sketch  the  proof. 

THEOREM  4.3.  Let  Assumptions  4.1  and  4.2  hold.  Then  if  Sq,  {r)n}>  and  h  are  sufficiently 
small,  and  {()„}  is  bounded  from  below,  then  xn  — >  x*  and  (4.1)  holds.  If  F  is  semismooth  of  order 
p,  then  (4.2)  holds. 

Proof.  Let 

n— 1 

tn  ^  '  Sn. 

n= 0 

The  SER  formula  implies  that 

Sn<S0/\\F(xn)\\. 

For  now  assume  that  F  is  semismooth  of  order  p  and  that  r/n  =  0  for  all  n.  Our  assumptions 
imply  that  there  is  T  such  that  x(t)  £  B(x*,  6^/2),  the  ball  of  local  convergence  from  Theorem  4.1, 
for  all  t  >  T.  Lemma  4.2  implies  that  for  small  S,  the  Ttc  iteration  is  an  accurate  integrator  for 

(2.8)  in  the  sense  that 

(4.9)  §£„  -  x(tn) ||  =  0(SP  +  nh), 

where  S  =  max0<k<n  Sk-  Hence,  we  can  select  S0  and  h  such  that  xn  £  S(eG)  until  tn  >  T. 

If  F  is  semismooth  and  {pn}  is  non-zero,  then  (4.9)  becomes  (see  [18]) 

n 

(4.10)  \\xn  -  x(tn) ||  =  0(nh  +  Y'fijSj)  +  o(l),  as  S  — >•  0  , 

j=o 

and  the  convergence  result  still  holds  if,  say  r/n  =  ()(do).  □ 

5.  Numerical  Example.  We  illustrate  the  results  with  a  simple  one-dimensional  example 
taken  from  [1,2,4].  This  example  is  sufficient  to  illustrate  the  convergence  results  in  this  paper, 
and  allows  us  to  refine  the  grids  to  a  degree  that  was  not  possible  in  the  two  and  three  dimensional 
results  that  motivated  this  paper  [10, 14, 15,24,30,31]. 

We  use  direct  methods  to  compute  the  Newton  step  in  this  section,  so  r/n  =  0.  In  all  but  §  5.2, 
we  compute  V  £  OF (x)  analytically,  so  h  =  0  in  those  computations. 

This  example,  taken  from  [4],  is  a  Lipschitz  reformulation  of  the  boundary  value  problem  [1,2] 


uZz  +  Amax(0,«)p  =  0,  z  £  (0, 1), 
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with  boundary  data 

■u(O)  =  u(  1)  =  0. 

and  p  e  (0, 1). 

The  reformulation  adds  a  new  variable 


v  = 


up  if  u  >  0 

u  if  u  <  0 


to  obtain  a  Lipschitz  continuous 

(5.1)  F(x)  =  ( 

where 


elliptic-algebraic  system,  F(x)  =  0, 


f{u,v) 

\  (  ~uzz  +  A  nrax(0,  v 

g(u,v) 

)  -  V  u-uj{v) 

/  \ 

(  v xlp  if  v  >  0 

UJ{V) 

|  v  if  v  <  0 

where  x 


(u,  v)T  and 


We  report  results  on  a  pseudo-transient  continuation  (\Ptc  )  approach.  We  use  the  method 
from  [7]  which  is  designed  for  differential  algebraic  equations  (DAEs). 

The  reason  we  formulate  the  problem  with  DAE  (rather  than  ODE)  dynamics  is  that  the 
pseudo-time  variable  should  not  be  added  to  both  equations  in  (5.1),  but  only  the  first.  The  reason 
for  this  is  that  the  true  time-dependent  system  is 


ut  =  uzz  —  A  nrax(0,  u)p . 


and  that  the  auxiliary  variable  v  is  used  only  to  make  the  nonlinearity  Lipschitz  continuous.  One 
might  think  that  an  ODE  formulation  would  work  equally  well,  but,  in  fact,  the  ODE  formulation, 
which  does  not  model  the  physics,  failed  to  converge  in  our  testing. 

We  discretize  the  problem  with  a  central  differences,  using  a  difference  increment  of  5Z.  The 
nonsmooth  nonlinearity  is  a  substitution  operator,  and  its  generalized  Jacobian  is  a  set  of  diagonal 
matrices. 

We  report  on  several  computations  with  p  —  .1  and  A  =  200.  This  choice  leads  to  a  large 
“dead  core”  [1,2],  a  region  in  which  the  solution  vanishes.  We  plot  the  solution  in  Figure  5.1. 


Eig.  5.1.  Solution 
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S(i  —  1  and  8max  =  106  for  all  the  computations.  We  terminate  the  nonlinear  iteration  when 
either 

(5.2)  ||F(xn)||/||F(x0)||  <  10~13  or  ||sn||  <  lO’10, 

where  sn  =  xn+  \  —  xn.  In  the  tables  we  see  the  superlinear  convergence  clearly  in  the  reduction 
in  the  norms  of  the  steps,  this  is  consistent  with  the  estimate  sn  =  —en  +  o(||en||)  which  follows 
from  local  superlinear  convergence.  The  superlinear  convergence  is  less  visible  in  the  residual 
norms,  because  the  generalized  Jacobians  become  more  ill-conditioned  as  the  mesh  is  refined.  The 
residual  norms  begin  to  stagnate  after  a  reduction  of  10 12 . 


5.1.  Exact  Compuation  of  the  Generalized  Jacobian.  For  the  results  in  this  section  we  com¬ 
pute  the  generalized  Jacobian  analytically.  If  we  let  L$z  be  the  discretized  Laplacian,  we  can  write 


(5.3) 


F(x) 


f{u,v) 
g(u,v ) 


~LSzu  \  + 

u  —  v  —  max( 0,  v1^)  J 


X 

1 


max(0,  v), 


and  use  the  known  result  for  the  scalar  function  max(0,  v) 


(  0,  if  v  <  0 
<9max(0,  u)  =  <  [0, 1],  if v  —  0 
[  1,  if  v  >  0, 

we  obtain 

(5'4)  SF=(“f‘  -l-(l/P)m<L(0.^-^))  +  (o  t)SmaX(°’1’)' 

The  calculations  in  this  section  used  V(xn )  G  dF(xn).  We  may  use  any  choice  from  the 
set-valued  map  <9(0,  v)  and  we  choose  W  G  OF  using 


0,  if  v  <  0 
1,  if  v  >  0, 


G  <9max(0,n). 


With  this  choice,  the  Wvv  is  nonsingular.  Had  we  used  1  when  v  —  0,  the  Wvv  would  be  singular 
at  v  =  0. 

In  Figure  5.2  we  plot  the  norms  of  the  steps  and  nonlinear  residuals  together  with  the  growth 
of  5  for  a  mesh  of  width  Sz  =  1/2048.  5  grows  smoothly  in  the  early  phase  of  the  iteration  and 
reaches  its  maximum  rapidly.  The  superlinear  convergence  is  clearly  visible  in  the  curve  for  the 
norms  of  the  steps.  The  Jacobian  of  the  nonlinear  residual  has  a  condition  number  of  0(1/ h2),  and 
reflects  the  error  less  accurately. 


5.2.  Compuation  of  the  Generalized  Jacobian  by  Differences.  For  the  results  in  this  section 
we  compute  the  generalized  Jacobian  with  several  choices  of  difference.  The  results  were  similar 
for  all  the  meshes.  In  Figures  5.2  and  5.3  we  plot  residual  and  step  norm  histories  for 

•  analytic  generalized  Jacobian  (Exact), 

•  forward  differences,  increment  10-8  (F-8),  and 

•  forward  differences,  increment  10  10  (F-10), 

for  a  mesh  of  width  Sz  =  1/2048  and  20  iterations.  In  this  way  we  can  clearly  see  the  point  at 
which  the  iteration  stagnates.  As  we  predicted  in  §  3.3,  the  iteration  is  more  accurate  when  the 
difference  increment  is  10~10  ~  e2//3  than  it  is  with  the  standard  choice  of  10-8  ~  e:1  2. 
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LlG.  5.2.  Analytic  Generalized  Jacobian 


LlG.  5.3.  Norms  of  the  Steps  and  Residuals. 


5.3.  Mesh  Dependence  and  Nested  Iteration.  We  used  the  analytic  d F  (5.4)  in  the  compu¬ 
tations  reported  in  this  section. 

In  Figure  5.4  we  plot  the  progress  of  the  iteration  for  mesh  sizes  of  1/128, 1/512,  and  1/2048, 
terminating  the  iteration  when  ||s||  <  10-13.  In  this  way  we  can  examine  the  dependence  of  the 
convergence  on  the  mesh  width.  While  the  convergence  in  the  early  phase  is  identical  for  the  three 
meshes,  and  superlinear  in  the  terminal  phase,  the  global  convergence  becomes  slower  as  the  mesh 
is  refined. 

Nested  iteration  or  grid  sequencing  means  to  solve  the  problem  to  high  precision  on  a  coarse 
mesh,  interpolate  to  a  finer  mesh  in  such  a  way  that  the  interpolation  error  can  be  corrected  with  a 
few  (eg  one)  iterations,  and  to  continue  this  until  one  has  a  solution  on  a  target,  finest  mesh.  We 
set  5  =  106  for  the  finer  meshes,  under  the  assumption  that  we  are  in  the  locally  convergent  phase 
of  the  iteration. 
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Fig.  5.4.  Mesh-dependence  of  Convergence 


Iterations 


For  this  example,  one  would  hope  not  only  to  eliminate  the  mesh-dependency  in  the  iteration 
history  that  are  visible  in  Figure  5.4,  but  also  to  approximate  the  solution  up  to  truncation  error  at 
each  level. 

This  was  a  successful  strategy.  However,  the  results  must  be  interpreted  in  light  of  the  conti¬ 
nuity  properties  of  the  solution,  u*  is  Lipschitz  continuously  differentiable  on  [0,1],  but  if  p  <  1/2, 
v*  =  (u*  )p  is  not.  This  means  that  linear  interpolation  will  not  approximate  v*  to  second  order  if 
p  <  1/2.  To  address  this  we  interpolate  u  from  the  coarse  to  fine  mesh  with  linear  interpolation, 
and  then  compute  v  as 

v  =  max(  0,  u)p . 

In  tables  5.1  and  5.2  we  report  the  residual  and  step  norms  on  a  sequence  of  meshes  {2_n}4=6l 
for  p  =  .1  and  p  =  .5.  The  initial  steps  at  each  mesh  reflect  both  the  error  in  the  initial  iterate  and 
the  truncation  error  in  the  interpolation. 

The  iterations  for  both  values  of  p  show  that  we  have  recovered  mesh  independence  in  the 
sense  that  the  iteration  requires  a  roughly  constant  number  of  steps  to  terminate  at  each  level.  The 
tables  for  p  =  1/2  clearly  show  second  order  convergence.  The  interpolation  error  for  p  —  .1  is 
visible  in  the  sizes  of  the  initial  steps. 


Table  5.1 

Step  norms:  p  =  .1,  nested  iteration 


n\5z 

1/64 

1/128 

1/256 

1/512 

1/1024 

1/2048 

0 

4.20e+00 

2.02e-02 

1.02e-02 

5.72e-03 

3.45e-03 

3.61e-03 

1 

3.53e+00 

1.13e-02 

1.23e-03 

1.13e-03 

2.14e-03 

6.16e-04 

2 

3.91e-02 

8.95e-04 

1.56e-04 

2.15e-04 

1.58e-04 

1.37e-05 

3 

4.11e-03 

6.44e-05 

2.19e-06 

6.18e-06 

7.24e-05 

7.14e-08 

4 

6.89e-04 

3.26e-07 

2.23e-12 

5 

1.94e-05 

6 

1.47e-08 
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Table  5.2 

Step  norms:  p  =  .5,  nested  iteration 


n\5z 

1/64 

1/128 

1/256 

1/512 

1/1024 

1/2048 

0 

1.32e+00 

1.52e-03 

3.87e-04 

9.74e-05 

2.44e-05 

6.13e-06 

1 

5.29e-01 

4.37e-05 

7.89e-06 

1.39e-06 

2.47e-07 

4.25e-08 

2 

5.20e-03 

1.10e-06 

9.73e-08 

4.75e-08 

4.21e-09 

3 

6.59e-05 

3.83e-09 

1.76e-ll 

4 

2.73e-05 

5 

9.74e-08 
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